!-------------------------------------------------------------------------------

! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2020 EDF S.A.
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 2 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.

!-------------------------------------------------------------------------------

!> \file chem_luscheme2.f90
!> \brief Routines for atmospheric chemical scheme 2
!> \remarks
!>  These routines are automatically generated by SPACK
!>  See CEREA: http://cerea.enpc.fr/polyphemus

!> kinetic_2
!> \brief Computation of kinetic rates for atmospheric chemistry
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name              role
!------------------------------------------------------------------------------
!> \param[in]     nr                total number of chemical reactions
!> \param[in]     option_photolysis flag to activate or not photolysis reactions
!> \param[in]     azi               solar zenith angle
!> \param[in]     att               atmospheric attenuation variable
!> \param[in]     temp              temperature
!> \param[in]     press             pressure
!> \param[in]     xlw               water massic fraction
!> \param[out]    rk(nr)            kinetic rates
!______________________________________________________________________________

subroutine kinetic_2(nr,rk,temp,xlw,press,azi,att,                  &
     option_photolysis)

implicit none

! Arguments

integer nr
double precision rk(nr),temp,xlw,press
double precision azi, att
integer option_photolysis

! Local variables

double precision effko,rapk,summ
double precision ylh2o

!     Compute third body.
!     Conversion = Avogadro*1d-6/Perfect gas constant.
!     PRESS in Pascal, SUMM in molecules/cm3, TEMP in Kelvin

summ = press * 7.243d16 / temp

!     Number of water molecules computed from the massic fraction
!     (absolute humidity)

ylh2o = 29.d0*summ*xlw/(18.d0+11.d0*xlw)

!     For the zenithal angle at tropics

azi=abs(azi)

rk(  1) =  dexp(-0.8860689615829534d+02                           &
 - ( -0.5300000000000000d+03 )/temp)
rk(  1) = rk(  1) * summ * 0.2d0
rk(  2) =  dexp(-0.2653240882726044d+02                           &
 - (  0.1500000000000000d+04 )/temp)

if(option_photolysis.eq.2) then
 rk(  3)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk(  3)=-0.1302720567168795d-07
 rk(  3)=-0.7822279432831311d-06+(azi- 0.00d+00) * rk(  3)
 rk(  3)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk(  3)
 rk(  3)= 0.9310260000000001d-02+(azi- 0.00d+00) * rk(  3)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk(  3)= 0.3771617015067078d-08
 rk(  3)=-0.1173044113433769d-05+(azi- 0.10d+02) * rk(  3)
 rk(  3)=-0.1955272056716901d-04+(azi- 0.10d+02) * rk(  3)
 rk(  3)= 0.9219010000000000d-02+(azi- 0.10d+02) * rk(  3)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk(  3)=-0.5859262388581815d-08
 rk(  3)=-0.1059895602981758d-05+(azi- 0.20d+02) * rk(  3)
 rk(  3)=-0.4188211773132428d-04+(azi- 0.20d+02) * rk(  3)
 rk(  3)= 0.8909950000000000d-02+(azi- 0.20d+02) * rk(  3)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk(  3)=-0.7024567460738029d-08
 rk(  3)=-0.1235673474639213d-05+(azi- 0.30d+02) * rk(  3)
 rk(  3)=-0.6483780850753392d-04+(azi- 0.30d+02) * rk(  3)
 rk(  3)= 0.8379279999999999d-02+(azi- 0.30d+02) * rk(  3)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk(  3)=-0.9202467768466835d-08
 rk(  3)=-0.1446410498461367d-05+(azi- 0.40d+02) * rk(  3)
 rk(  3)=-0.9165864823853972d-04+(azi- 0.40d+02) * rk(  3)
 rk(  3)= 0.7600310000000000d-02+(azi- 0.40d+02) * rk(  3)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk(  3)=-0.1612556146540100d-07
 rk(  3)=-0.1722484531515342d-05+(azi- 0.50d+02) * rk(  3)
 rk(  3)=-0.1233475985383066d-03+(azi- 0.50d+02) * rk(  3)
 rk(  3)= 0.6529880000000000d-02+(azi- 0.50d+02) * rk(  3)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk(  3)= 0.3226471363007382d-07
 rk(  3)=-0.2206251375477548d-05+(azi- 0.60d+02) * rk(  3)
 rk(  3)=-0.1626349576082332d-03+(azi- 0.60d+02) * rk(  3)
 rk(  3)= 0.5108030000000000d-02+(azi- 0.60d+02) * rk(  3)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk(  3)= 0.2027078243961372d-06
 rk(  3)=-0.1238309966574737d-05+(azi- 0.70d+02) * rk(  3)
 rk(  3)=-0.1970805710287543d-03+(azi- 0.70d+02) * rk(  3)
 rk(  3)= 0.3293320000000000d-02+(azi- 0.70d+02) * rk(  3)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk(  3)=-0.7448311471194499d-07
 rk(  3)= 0.3626677818932555d-05+(azi- 0.78d+02) * rk(  3)
 rk(  3)=-0.1779736282099126d-03+(azi- 0.78d+02) * rk(  3)
 rk(  3)= 0.1741210000000000d-02+(azi- 0.78d+02) * rk(  3)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk(  3)= 0.2490309929270573d-05
 rk(  3)= 0.1839083065842406d-05+(azi- 0.86d+02) * rk(  3)
 rk(  3)=-0.1342475411316713d-03+(azi- 0.86d+02) * rk(  3)
 rk(  3)= 0.5113930000000000d-03+(azi- 0.86d+02) * rk(  3)
elseif(azi.ge.90.d0)then
 rk(  3)= 0.1632080000000000d-03
endif
if(att.lt.0.99999) rk(  3) = rk(  3) * att
endif

rk(  4) = summ * 6.0d-34 * (temp/3.d2) ** (-2.4d0)
rk(  4) = rk(  4) * summ * 0.2d0

if(option_photolysis.eq.2) then
 rk(  5)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk(  5)=-0.5928286648807996d-09
 rk(  5)=-0.3096171335119280d-07+(azi- 0.00d+00) * rk(  5)
 rk(  5)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk(  5)
 rk(  5)= 0.4927580000000000d-03+(azi- 0.00d+00) * rk(  5)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk(  5)= 0.1444859946426876d-09
 rk(  5)=-0.4874657329761677d-07+(azi- 0.10d+02) * rk(  5)
 rk(  5)=-0.7970828664880956d-06+(azi- 0.10d+02) * rk(  5)
 rk(  5)= 0.4890690000000000d-03+(azi- 0.10d+02) * rk(  5)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk(  5)=-0.5511531369012520d-10
 rk(  5)=-0.4441199345833616d-07+(azi- 0.20d+02) * rk(  5)
 rk(  5)=-0.1728668534047625d-05+(azi- 0.20d+02) * rk(  5)
 rk(  5)= 0.4763680000000000d-03+(azi- 0.20d+02) * rk(  5)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk(  5)=-0.3000247398821449d-09
 rk(  5)=-0.4606545286904014d-07+(azi- 0.30d+02) * rk(  5)
 rk(  5)=-0.2633442997321385d-05+(azi- 0.30d+02) * rk(  5)
 rk(  5)= 0.4545850000000000d-03+(azi- 0.30d+02) * rk(  5)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk(  5)=-0.2397857267812366d-09
 rk(  5)=-0.5506619506550444d-07+(azi- 0.40d+02) * rk(  5)
 rk(  5)=-0.3644759476666826d-05+(azi- 0.40d+02) * rk(  5)
 rk(  5)= 0.4233440000000000d-03+(azi- 0.40d+02) * rk(  5)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk(  5)=-0.1844832352993067d-08
 rk(  5)=-0.6225976686893995d-07+(azi- 0.50d+02) * rk(  5)
 rk(  5)=-0.4818019096011278d-05+(azi- 0.50d+02) * rk(  5)
 rk(  5)= 0.3811500000000000d-03+(azi- 0.50d+02) * rk(  5)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk(  5)= 0.1101151387530619d-09
 rk(  5)=-0.1176047374587370d-06+(azi- 0.60d+02) * rk(  5)
 rk(  5)=-0.6616664139287950d-05+(azi- 0.60d+02) * rk(  5)
 rk(  5)= 0.3248990000000000d-03+(azi- 0.60d+02) * rk(  5)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk(  5)=-0.1557211541866474d-07
 rk(  5)=-0.1143012832961418d-06+(azi- 0.70d+02) * rk(  5)
 rk(  5)=-0.8935724346837023d-05+(azi- 0.70d+02) * rk(  5)
 rk(  5)= 0.2470820000000000d-03+(azi- 0.70d+02) * rk(  5)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk(  5)= 0.4048472604232427d-07
 rk(  5)=-0.4880320533439059d-06+(azi- 0.78d+02) * rk(  5)
 rk(  5)=-0.1375439103995816d-04+(azi- 0.78d+02) * rk(  5)
 rk(  5)= 0.1603080000000000d-03+(azi- 0.78d+02) * rk(  5)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk(  5)= 0.2066880316654862d-06
 rk(  5)= 0.4836013716715513d-06+(azi- 0.86d+02) * rk(  5)
 rk(  5)=-0.1378983649333310d-04+(azi- 0.86d+02) * rk(  5)
 rk(  5)= 0.3976700000000000d-04+(azi- 0.86d+02) * rk(  5)
elseif(azi.ge.90.d0)then
 rk(  5)= 0.5573310000000000d-05
endif
if(att.lt.0.99999) rk(  5) = rk(  5) * att
endif


if(option_photolysis.eq.2) then
 rk(  6)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk(  6)=-0.8776629099833464d-10
 rk(  6)=-0.1165033709001661d-07+(azi- 0.00d+00) * rk(  6)
 rk(  6)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk(  6)
 rk(  6)= 0.3523480000000000d-04+(azi- 0.00d+00) * rk(  6)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk(  6)= 0.1474988729949909d-09
 rk(  6)=-0.1428332581996665d-07+(azi- 0.10d+02) * rk(  6)
 rk(  6)=-0.2593366290998327d-06+(azi- 0.10d+02) * rk(  6)
 rk(  6)= 0.3398200000000000d-04+(azi- 0.10d+02) * rk(  6)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk(  6)= 0.1300707990183827d-09
 rk(  6)=-0.9858359630116941d-08+(azi- 0.20d+02) * rk(  6)
 rk(  6)=-0.5007534836006686d-06+(azi- 0.20d+02) * rk(  6)
 rk(  6)= 0.3010780000000000d-04+(azi- 0.20d+02) * rk(  6)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk(  6)= 0.1988179309314732d-09
 rk(  6)=-0.5956235659565481d-08+(azi- 0.30d+02) * rk(  6)
 rk(  6)=-0.6588994364974921d-06+(azi- 0.30d+02) * rk(  6)
 rk(  6)= 0.2424450000000000d-04+(azi- 0.30d+02) * rk(  6)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk(  6)= 0.2219574772557277d-09
 rk(  6)= 0.8302268378835231d-11+(azi- 0.40d+02) * rk(  6)
 rk(  6)=-0.7183787704093613d-06+(azi- 0.40d+02) * rk(  6)
 rk(  6)= 0.1725870000000000d-04+(azi- 0.40d+02) * rk(  6)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk(  6)= 0.1913521600455895d-09
 rk(  6)= 0.6667026586050136d-08+(azi- 0.50d+02) * rk(  6)
 rk(  6)=-0.6516254818650674d-06+(azi- 0.50d+02) * rk(  6)
 rk(  6)= 0.1029770000000000d-04+(azi- 0.50d+02) * rk(  6)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk(  6)= 0.1602388256152816d-10
 rk(  6)= 0.1240759138741867d-07+(azi- 0.60d+02) * rk(  6)
 rk(  6)=-0.4608793021303539d-06+(azi- 0.60d+02) * rk(  6)
 rk(  6)= 0.4639500000000000d-05+(azi- 0.60d+02) * rk(  6)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk(  6)=-0.3089359890783960d-09
 rk(  6)= 0.1288830786428400d-07+(azi- 0.70d+02) * rk(  6)
 rk(  6)=-0.2079203096133002d-06+(azi- 0.70d+02) * rk(  6)
 rk(  6)= 0.1287490000000000d-05+(azi- 0.70d+02) * rk(  6)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk(  6)=-0.2034952628632162d-09
 rk(  6)= 0.5473844126395715d-08+(azi- 0.78d+02) * rk(  6)
 rk(  6)=-0.6102309368797090d-07+(azi- 0.78d+02) * rk(  6)
 rk(  6)= 0.2908040000000000d-06+(azi- 0.78d+02) * rk(  6)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk(  6)= 0.1623544916128788d-09
 rk(  6)= 0.5899578175158973d-09+(azi- 0.86d+02) * rk(  6)
 rk(  6)=-0.1251267813581064d-07+(azi- 0.86d+02) * rk(  6)
 rk(  6)= 0.4875570000000000d-07+(azi- 0.86d+02) * rk(  6)
elseif(azi.ge.90.d0)then
 rk(  6)= 0.1853500000000000d-07
endif
if(att.lt.0.99999) rk(  6) = rk(  6) * att
endif

rk(  7) =  dexp(-0.2458649867820512d+02                           &
 - ( -0.1020000000000000d+03 )/temp)
rk(  7) = rk(  7) * summ
rk(  8) =  0.2200000000000000d-09
rk(  8) = rk(  8) * ylh2o

if(option_photolysis.eq.2) then
 rk(  9)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk(  9)=-0.2887225450832658d-08
 rk(  9)=-0.1810277454916759d-06+(azi- 0.00d+00) * rk(  9)
 rk(  9)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk(  9)
 rk(  9)= 0.2046710000000000d-02+(azi- 0.00d+00) * rk(  9)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk(  9)= 0.8216763524985595d-09
 rk(  9)=-0.2676445090166556d-06+(azi- 0.10d+02) * rk(  9)
 rk(  9)=-0.4486722545083316d-05+(azi- 0.10d+02) * rk(  9)
 rk(  9)= 0.2025720000000000d-02+(azi- 0.10d+02) * rk(  9)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk(  9)=-0.1309479959161308d-08
 rk(  9)=-0.2429942184416991d-06+(azi- 0.20d+02) * rk(  9)
 rk(  9)=-0.9593109819666860d-05+(azi- 0.20d+02) * rk(  9)
 rk(  9)= 0.1954910000000000d-02+(azi- 0.20d+02) * rk(  9)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk(  9)=-0.1523756515853649d-08
 rk(  9)=-0.2822786172165394d-06+(azi- 0.30d+02) * rk(  9)
 rk(  9)=-0.1484583817624924d-04+(azi- 0.30d+02) * rk(  9)
 rk(  9)= 0.1833370000000000d-02+(azi- 0.30d+02) * rk(  9)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk(  9)=-0.1745493977424016d-08
 rk(  9)=-0.3279913126921461d-06+(azi- 0.40d+02) * rk(  9)
 rk(  9)=-0.2094853747533609d-04+(azi- 0.40d+02) * rk(  9)
 rk(  9)= 0.1655160000000000d-02+(azi- 0.40d+02) * rk(  9)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk(  9)=-0.2754267574450560d-08
 rk(  9)=-0.3803561320148624d-06+(azi- 0.50d+02) * rk(  9)
 rk(  9)=-0.2803201192240625d-04+(azi- 0.50d+02) * rk(  9)
 rk(  9)= 0.1411130000000000d-02+(azi- 0.50d+02) * rk(  9)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk(  9)= 0.1037656427523324d-07
 rk(  9)=-0.4629841592484096d-06+(azi- 0.60d+02) * rk(  9)
 rk(  9)=-0.3646541483503878d-04+(azi- 0.60d+02) * rk(  9)
 rk(  9)= 0.1090020000000000d-02+(azi- 0.60d+02) * rk(  9)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk(  9)= 0.4335158727139456d-07
 rk(  9)=-0.1516872309913989d-06+(azi- 0.70d+02) * rk(  9)
 rk(  9)=-0.4261212873743828d-04+(azi- 0.70d+02) * rk(  9)
 rk(  9)= 0.6894440000000000d-03+(azi- 0.70d+02) * rk(  9)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk(  9)=-0.2976581610911796d-07
 rk(  9)= 0.8887508635220705d-06+(azi- 0.78d+02) * rk(  9)
 rk(  9)=-0.3671561967719464d-04+(azi- 0.78d+02) * rk(  9)
 rk(  9)= 0.3610350000000000d-03+(azi- 0.78d+02) * rk(  9)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk(  9)= 0.5586598403857848d-06
 rk(  9)= 0.1743712769036732d-06+(azi- 0.86d+02) * rk(  9)
 rk(  9)=-0.2821064255377481d-04+(azi- 0.86d+02) * rk(  9)
 rk(  9)= 0.1089500000000000d-03+(azi- 0.86d+02) * rk(  9)
elseif(azi.ge.90.d0)then
 rk(  9)= 0.3465160000000000d-04
endif
if(att.lt.0.99999) rk(  9) = rk(  9) * att
endif


if(option_photolysis.eq.2) then
 rk( 10)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk( 10)=-0.1441479345432039d-10
 rk( 10)=-0.1242452065456794d-08+(azi- 0.00d+00) * rk( 10)
 rk( 10)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 10)
 rk( 10)= 0.8394580000000000d-05+(azi- 0.00d+00) * rk( 10)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk( 10)= 0.8244380362960478d-11
 rk( 10)=-0.1674895869086405d-08+(azi- 0.10d+02) * rk( 10)
 rk( 10)=-0.2917347934543199d-07+(azi- 0.10d+02) * rk( 10)
 rk( 10)= 0.8255920000000000d-05+(azi- 0.10d+02) * rk( 10)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk( 10)= 0.1172720024787734d-12
 rk( 10)=-0.1427564458197592d-08+(azi- 0.20d+02) * rk( 10)
 rk( 10)=-0.6019808261827194d-07+(azi- 0.20d+02) * rk( 10)
 rk( 10)= 0.7804940000000000d-05+(azi- 0.20d+02) * rk( 10)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk( 10)= 0.4506531627127201d-11
 rk( 10)=-0.1424046298123240d-08+(azi- 0.30d+02) * rk( 10)
 rk( 10)=-0.8871419018148013d-07+(azi- 0.30d+02) * rk( 10)
 rk( 10)= 0.7060320000000000d-05+(azi- 0.30d+02) * rk( 10)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk( 10)= 0.1086660148900755d-10
 rk( 10)=-0.1288850349309390d-08+(azi- 0.40d+02) * rk( 10)
 rk( 10)=-0.1158431566558070d-06+(azi- 0.40d+02) * rk( 10)
 rk( 10)= 0.6035280000000000d-05+(azi- 0.40d+02) * rk( 10)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk( 10)= 0.1803706241686108d-10
 rk( 10)=-0.9628523046392104d-09+(azi- 0.50d+02) * rk( 10)
 rk( 10)=-0.1383601831952927d-06+(azi- 0.50d+02) * rk( 10)
 rk( 10)= 0.4758830000000000d-05+(azi- 0.50d+02) * rk( 10)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk( 10)= 0.7329514884355585d-10
 rk( 10)=-0.4217404321336693d-09+(azi- 0.60d+02) * rk( 10)
 rk( 10)=-0.1522061105630206d-06+(azi- 0.60d+02) * rk( 10)
 rk( 10)= 0.3296980000000000d-05+(azi- 0.60d+02) * rk( 10)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk( 10)= 0.6172011386267615d-10
 rk( 10)= 0.1777114033174383d-08+(azi- 0.70d+02) * rk( 10)
 rk( 10)=-0.1386523745526241d-06+(azi- 0.70d+02) * rk( 10)
 rk( 10)= 0.1806040000000000d-05+(azi- 0.70d+02) * rk( 10)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk( 10)=-0.8279704635596216d-10
 rk( 10)= 0.3258396765877763d-08+(azi- 0.78d+02) * rk( 10)
 rk( 10)=-0.9836828816022045d-07+(azi- 0.78d+02) * rk( 10)
 rk( 10)= 0.8421570000000000d-06+(azi- 0.78d+02) * rk( 10)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk( 10)= 0.1082517324572734d-08
 rk( 10)= 0.1271267653334671d-08+(azi- 0.86d+02) * rk( 10)
 rk( 10)=-0.6213097280649386d-07+(azi- 0.86d+02) * rk( 10)
 rk( 10)= 0.2213560000000000d-06+(azi- 0.86d+02) * rk( 10)
elseif(azi.ge.90.d0)then
 rk( 10)= 0.6245350000000000d-07
endif
if(att.lt.0.99999) rk( 10) = rk( 10) * att
endif


if(option_photolysis.eq.2) then
 rk( 11)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk( 11)=-0.5816837387219519d-10
 rk( 11)=-0.5184316261278087d-08+(azi- 0.00d+00) * rk( 11)
 rk( 11)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 11)
 rk( 11)= 0.3220560000000000d-04+(azi- 0.00d+00) * rk( 11)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk( 11)= 0.3450512161659949d-10
 rk( 11)=-0.6929367477443941d-08+(azi- 0.10d+02) * rk( 11)
 rk( 11)=-0.1211368373872203d-06+(azi- 0.10d+02) * rk( 11)
 rk( 11)= 0.3162900000000000d-04+(azi- 0.10d+02) * rk( 11)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk( 11)= 0.3647887405787883d-11
 rk( 11)=-0.5894213828945960d-08+(azi- 0.20d+02) * rk( 11)
 rk( 11)=-0.2493726504511192d-06+(azi- 0.20d+02) * rk( 11)
 rk( 11)= 0.2975920000000000d-04+(azi- 0.20d+02) * rk( 11)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk( 11)= 0.2530332876025029d-10
 rk( 11)=-0.5784777206772343d-08+(azi- 0.30d+02) * rk( 11)
 rk( 11)=-0.3661625608083018d-06+(azi- 0.30d+02) * rk( 11)
 rk( 11)= 0.2667970000000000d-04+(azi- 0.30d+02) * rk( 11)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk( 11)= 0.6013879755321277d-10
 rk( 11)=-0.5025677343964861d-08+(azi- 0.40d+02) * rk( 11)
 rk( 11)=-0.4742671063156733d-06+(azi- 0.40d+02) * rk( 11)
 rk( 11)= 0.2246490000000000d-04+(azi- 0.40d+02) * rk( 11)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk( 11)= 0.1004414810269593d-09
 rk( 11)=-0.3221513417368319d-08+(azi- 0.50d+02) * rk( 11)
 rk( 11)=-0.5567390139290089d-06+(azi- 0.50d+02) * rk( 11)
 rk( 11)= 0.1727980000000000d-04+(azi- 0.50d+02) * rk( 11)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk( 11)= 0.3107352783389442d-09
 rk( 11)=-0.2082689865616575d-09+(azi- 0.60d+02) * rk( 11)
 rk( 11)=-0.5910368379682989d-06+(azi- 0.60d+02) * rk( 11)
 rk( 11)= 0.1149070000000000d-04+(azi- 0.60d+02) * rk( 11)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk( 11)= 0.1345128013872392d-09
 rk( 11)= 0.9113789363612173d-08+(azi- 0.70d+02) * rk( 11)
 rk( 11)=-0.5019816341977565d-06+(azi- 0.70d+02) * rk( 11)
 rk( 11)= 0.5870240000000000d-05+(azi- 0.70d+02) * rk( 11)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk( 11)=-0.1960445509611151d-09
 rk( 11)= 0.1234209659691269d-07+(azi- 0.78d+02) * rk( 11)
 rk( 11)=-0.3303345465135304d-06+(azi- 0.78d+02) * rk( 11)
 rk( 11)= 0.2506540000000000d-05+(azi- 0.78d+02) * rk( 11)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk( 11)= 0.2279277828298341d-08
 rk( 11)= 0.7637027373778166d-08+(azi- 0.86d+02) * rk( 11)
 rk( 11)=-0.1705015547476783d-06+(azi- 0.86d+02) * rk( 11)
 rk( 11)= 0.5533830000000000d-06+(azi- 0.86d+02) * rk( 11)
elseif(azi.ge.90.d0)then
 rk( 11)= 0.1394430000000000d-06
endif
if(att.lt.0.99999) rk( 11) = rk( 11) * att
endif


if(option_photolysis.eq.2) then
 rk( 12)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk( 12)=-0.7261045436751545d-10
 rk( 12)=-0.5576895456324842d-08+(azi- 0.00d+00) * rk( 12)
 rk( 12)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 12)
 rk( 12)= 0.4669460000000000d-04+(azi- 0.00d+00) * rk( 12)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk( 12)= 0.2853136310254372d-10
 rk( 12)=-0.7755209087350302d-08+(azi- 0.10d+02) * rk( 12)
 rk( 12)=-0.1333210454367515d-06+(azi- 0.10d+02) * rk( 12)
 rk( 12)= 0.4606430000000000d-04+(azi- 0.10d+02) * rk( 12)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk( 12)=-0.1901499804265726d-10
 rk( 12)=-0.6899268194273995d-08+(azi- 0.20d+02) * rk( 12)
 rk( 12)=-0.2798658182529944d-06+(azi- 0.20d+02) * rk( 12)
 rk( 12)= 0.4398410000000000d-04+(azi- 0.20d+02) * rk( 12)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk( 12)=-0.9471370931908783d-11
 rk( 12)=-0.7469718135553710d-08+(azi- 0.30d+02) * rk( 12)
 rk( 12)=-0.4235556815512711d-06+(azi- 0.30d+02) * rk( 12)
 rk( 12)= 0.4047650000000000d-04+(azi- 0.30d+02) * rk( 12)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk( 12)= 0.1440048177028887d-10
 rk( 12)=-0.7753859263510966d-08+(azi- 0.40d+02) * rk( 12)
 rk( 12)=-0.5757914555419174d-06+(azi- 0.40d+02) * rk( 12)
 rk( 12)= 0.3548450000000000d-04+(azi- 0.40d+02) * rk( 12)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk( 12)= 0.3856944385081257d-10
 rk( 12)=-0.7321844810402538d-08+(azi- 0.50d+02) * rk( 12)
 rk( 12)=-0.7265484962810543d-06+(azi- 0.50d+02) * rk( 12)
 rk( 12)= 0.2896560000000000d-04+(azi- 0.50d+02) * rk( 12)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk( 12)= 0.4136217428262900d-09
 rk( 12)=-0.6164761494878585d-08+(azi- 0.60d+02) * rk( 12)
 rk( 12)=-0.8614145593338596d-06+(azi- 0.60d+02) * rk( 12)
 rk( 12)= 0.2100650000000000d-04+(azi- 0.60d+02) * rk( 12)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk( 12)= 0.5376701572532502d-09
 rk( 12)= 0.6243890789914774d-08+(azi- 0.70d+02) * rk( 12)
 rk( 12)=-0.8606232663835604d-06+(azi- 0.70d+02) * rk( 12)
 rk( 12)= 0.1218950000000000d-04+(azi- 0.70d+02) * rk( 12)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk( 12)=-0.5508273899935273d-09
 rk( 12)= 0.1914797456399955d-07+(azi- 0.78d+02) * rk( 12)
 rk( 12)=-0.6574883435524898d-06+(azi- 0.78d+02) * rk( 12)
 rk( 12)= 0.5979410000000000d-05+(azi- 0.78d+02) * rk( 12)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk( 12)= 0.8530305661922505d-08
 rk( 12)= 0.5928117204100688d-08+(azi- 0.86d+02) * rk( 12)
 rk( 12)=-0.4568796094072541d-06+(azi- 0.86d+02) * rk( 12)
 rk( 12)= 0.1662950000000000d-05+(azi- 0.86d+02) * rk( 12)
elseif(azi.ge.90.d0)then
 rk( 12)= 0.4762210000000000d-06
endif
if(att.lt.0.99999) rk( 12) = rk( 12) * att
endif


if(option_photolysis.eq.2) then
 rk( 13)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk( 13)=-0.1292568102250478d-10
 rk( 13)=-0.1349143189774952d-08+(azi- 0.00d+00) * rk( 13)
 rk( 13)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 13)
 rk( 13)= 0.6105070000000000d-05+(azi- 0.00d+00) * rk( 13)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk( 13)= 0.1221704306751377d-10
 rk( 13)=-0.1736913620450095d-08+(azi- 0.10d+02) * rk( 13)
 rk( 13)=-0.3086056810225046d-07+(azi- 0.10d+02) * rk( 13)
 rk( 13)= 0.5957230000000000d-05+(azi- 0.10d+02) * rk( 13)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk( 13)= 0.7227508752451222d-11
 rk( 13)=-0.1370402328424685d-08+(azi- 0.20d+02) * rk( 13)
 rk( 13)=-0.6193372759099823d-07+(azi- 0.20d+02) * rk( 13)
 rk( 13)= 0.5487150000000000d-05+(azi- 0.20d+02) * rk( 13)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk( 13)= 0.1521292192268060d-10
 rk( 13)=-0.1153577065851153d-08+(azi- 0.30d+02) * rk( 13)
 rk( 13)=-0.8717352153375648d-07+(azi- 0.30d+02) * rk( 13)
 rk( 13)= 0.4738000000000000d-05+(azi- 0.30d+02) * rk( 13)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk( 13)= 0.2301080355683006d-10
 rk( 13)=-0.6971894081707449d-09+(azi- 0.40d+02) * rk( 13)
 rk( 13)=-0.1056811862739754d-06+(azi- 0.40d+02) * rk( 13)
 rk( 13)= 0.3766120000000000d-05+(azi- 0.40d+02) * rk( 13)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk( 13)= 0.3106386384999131d-10
 rk( 13)=-0.6865301465823343d-11+(azi- 0.50d+02) * rk( 13)
 rk( 13)=-0.1127217333703412d-06+(azi- 0.50d+02) * rk( 13)
 rk( 13)= 0.2662600000000000d-05+(azi- 0.50d+02) * rk( 13)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk( 13)= 0.4531874104319860d-10
 rk( 13)= 0.9250506140339690d-09+(azi- 0.60d+02) * rk( 13)
 rk( 13)=-0.1035398802446585d-06+(azi- 0.60d+02) * rk( 13)
 rk( 13)= 0.1565760000000000d-05+(azi- 0.60d+02) * rk( 13)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk( 13)=-0.1303761111962380d-10
 rk( 13)= 0.2284612845330880d-08+(azi- 0.70d+02) * rk( 13)
 rk( 13)=-0.7144324565100488d-07+(azi- 0.70d+02) * rk( 13)
 rk( 13)= 0.6681850000000000d-06+(azi- 0.70d+02) * rk( 13)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk( 13)=-0.4077238229824927d-10
 rk( 13)= 0.1971710178462450d-08+(azi- 0.78d+02) * rk( 13)
 rk( 13)=-0.3739266146067857d-07+(azi- 0.78d+02) * rk( 13)
 rk( 13)= 0.2361790000000000d-06+(azi- 0.78d+02) * rk( 13)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk( 13)= 0.1193377495821847d-09
 rk( 13)= 0.9931730033112436d-09+(azi- 0.86d+02) * rk( 13)
 rk( 13)=-0.1367359600643481d-07+(azi- 0.86d+02) * rk( 13)
 rk( 13)= 0.4235170000000000d-07+(azi- 0.86d+02) * rk( 13)
elseif(azi.ge.90.d0)then
 rk( 13)= 0.1118570000000000d-07
endif
if(att.lt.0.99999) rk( 13) = rk( 13) * att
endif

effko =  0.4900000000000000d-02* (temp / 3.d2)                    &
           **(- ( 0.0000000000000000d+00))*                       &
           dexp(- 0.1210000000000000d+05/temp)
rapk =  0.5400000000000000d+17* (temp / 3.d2)                     &
            **(- ( 0.0000000000000000d+00))*                      &
           dexp(- 0.1383000000000000d+05/temp)
rk( 14) = (effko*summ / ( 1.0d0 + effko * summ /                  &
         rapk)) * 0.3000d+00** (1.0d0 / (1.0d0 +                  &
          (dlog10(effko * summ / rapk))**2))
effko =  0.2700000000000000d-27* (temp / 3.d2)                    &
           **(- ( 0.7100000000000000d+01))
rapk =  0.1200000000000000d-10* (temp / 3.d2)                     &
            **(- ( 0.9000000000000000d+00))
rk( 15) = (effko * summ /                                         &
            ( 1.0d0 + effko * summ / rapk)) *                     &
            0.3000d+00** (1.0d0 / (1.0d0 +                        &
             (dlog10(effko * summ / rapk))**2))
rk( 16) =  dexp(-0.2553915705425015d+02                           &
 - ( -0.2700000000000000d+03 )/temp)
rk( 17) =  dexp(-0.2637825814743318d+02                           &
 - ( -0.2500000000000000d+03 )/temp)
effko =  0.7000000000000000d-30* (temp / 3.d2)                    &
           **(- ( 0.2600000000000000d+01))
rapk =  0.3600000000000000d-10* (temp / 3.d2)                     &
            **(- ( 0.1000000000000000d+00))
rk( 18) = (effko * summ /                                         &
            ( 1.0d0 + effko * summ / rapk)) *                     &
            0.6000d+00** (1.0d0 / (1.0d0 +                        &
             (dlog10(effko * summ / rapk))**2))
rk( 19) =  dexp(-0.8322449114623726d+01                           &
 + (-0.3000000000000000d+01 * dlog(temp)) )
effko =  0.3300000000000000d-30* (temp / 3.d2)                    &
           **(- ( 0.3300000000000000d+01))
rapk =  0.1500000000000000d-11* (temp / 3.d2)                     &
            **(- ( 0.0000000000000000d+00))
rk( 20) = (effko * summ /                                         &
            ( 1.0d0 + effko * summ / rapk)) *                     &
            0.6000d+00** (1.0d0 / (1.0d0 +                        &
             (dlog10(effko * summ / rapk))**2))
rk( 21) =  dexp(-0.2975128465212864d+02                           &
 - (  0.2450000000000000d+04 )/temp)
rk( 22) =  dexp(-0.2492297091482634d+02                           &
 - ( -0.1700000000000000d+03 )/temp)
rk( 23) =  dexp(-0.3073211390514037d+02                           &
 - (  0.1260000000000000d+04 )/temp)
effko =  0.2000000000000000d-29* (temp / 3.d2)                    &
           **(- ( 0.4400000000000000d+01))
rapk =  0.1400000000000000d-11* (temp / 3.d2)                     &
            **(- ( 0.7000000000000000d+00))
rk( 24) = (effko * summ /                                         &
            ( 1.0d0 + effko * summ / rapk)) *                     &
            0.6000d+00** (1.0d0 / (1.0d0 +                        &
             (dlog10(effko * summ / rapk))**2))
effko =  0.1300000000000000d-02* (temp / 3.d2)                    &
           **(- ( 0.3500000000000000d+01))*                       &
           dexp(- 0.1100000000000000d+05/temp)
rapk =  0.9700000000000000d+15* (temp / 3.d2)                     &
            **(- (-0.1000000000000000d+00))*                      &
           dexp(- 0.1108000000000000d+05/temp)
rk( 25) = (effko*summ / ( 1.0d0 + effko * summ /                  &
         rapk)) * 0.4500d+00** (1.0d0 / (1.0d0 +                  &
          (dlog10(effko * summ / rapk))**2))
rk( 26) =  0.1000000000000000d-21
rk( 26) = rk( 26) * ylh2o
rk( 27) =  dexp(-0.2667550967090111d+02                           &
 - ( -0.3650000000000000d+03 )/temp)
rk( 28) =  0.6800000000000000d-13
rk( 29) = 2.3d-13 * dexp(600.0d0 / temp)                          &
            + 1.7d-33* summ * dexp(1000.0d0 / temp)
rk( 30) = 3.22d-34 * dexp(2800.0d0 / temp) +                      &
            2.38d-54 * summ * dexp(3200.0d0 / temp)
rk( 30) = rk( 30) * ylh2o

if(option_photolysis.eq.2) then
 rk( 31)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk( 31)=-0.2293082537254139d-06
 rk( 31)=-0.1232691746274577d-04+(azi- 0.00d+00) * rk( 31)
 rk( 31)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 31)
 rk( 31)= 0.2395610000000000d+00+(azi- 0.00d+00) * rk( 31)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk( 31)=-0.2107523882380334d-07
 rk( 31)=-0.1920616507450818d-04+(azi- 0.10d+02) * rk( 31)
 rk( 31)=-0.3153308253725395d-03+(azi- 0.10d+02) * rk( 31)
 rk( 31)= 0.2380990000000000d+00+(azi- 0.10d+02) * rk( 31)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk( 31)= 0.2060920902067111d-07
 rk( 31)=-0.1983842223922230d-04+(azi- 0.20d+02) * rk( 31)
 rk( 31)=-0.7057766985098441d-03+(azi- 0.20d+02) * rk( 31)
 rk( 31)= 0.2330040000000000d+00+(azi- 0.20d+02) * rk( 31)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk( 31)=-0.2323615972588770d-06
 rk( 31)=-0.1922014596860219d-04+(azi- 0.30d+02) * rk( 31)
 rk( 31)=-0.1096362380588088d-02+(azi- 0.30d+02) * rk( 31)
 rk( 31)= 0.2239830000000000d+00+(azi- 0.30d+02) * rk( 31)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk( 31)=-0.2281628199851812d-06
 rk( 31)=-0.2619099388636854d-04+(azi- 0.40d+02) * rk( 31)
 rk( 31)=-0.1550473779137796d-02+(azi- 0.40d+02) * rk( 31)
 rk( 31)= 0.2108650000000000d+00+(azi- 0.40d+02) * rk( 31)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk( 31)=-0.1187987122800317d-05
 rk( 31)=-0.3303587848592447d-04+(azi- 0.50d+02) * rk( 31)
 rk( 31)=-0.2142742502860728d-02+(azi- 0.50d+02) * rk( 31)
 rk( 31)= 0.1925130000000000d+00+(azi- 0.50d+02) * rk( 31)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk( 31)=-0.3438886888129482d-06
 rk( 31)=-0.6867549216993743d-04+(azi- 0.60d+02) * rk( 31)
 rk( 31)=-0.3159856209419318d-02+(azi- 0.60d+02) * rk( 31)
 rk( 31)= 0.1665940000000000d+00+(azi- 0.60d+02) * rk( 31)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk( 31)=-0.7936212779117297d-05
 rk( 31)=-0.7899215283432154d-04+(azi- 0.70d+02) * rk( 31)
 rk( 31)=-0.4636532659461956d-02+(azi- 0.70d+02) * rk( 31)
 rk( 31)= 0.1277840000000000d+00+(azi- 0.70d+02) * rk( 31)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk( 31)= 0.2654257866666759d-04
 rk( 31)=-0.2694612595331436d-03+(azi- 0.78d+02) * rk( 31)
 rk( 31)=-0.7424159958401733d-02+(azi- 0.78d+02) * rk( 31)
 rk( 31)= 0.8157290000000000d-01+(azi- 0.78d+02) * rk( 31)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk( 31)= 0.7705999956653109d-04
 rk( 31)= 0.3675606284668786d-03+(azi- 0.86d+02) * rk( 31)
 rk( 31)=-0.6639365006930298d-02+(azi- 0.86d+02) * rk( 31)
 rk( 31)= 0.1852390000000000d-01+(azi- 0.86d+02) * rk( 31)
elseif(azi.ge.90.d0)then
 rk( 31)= 0.2779250000000000d-02
endif
if(att.lt.0.99999) rk( 31) = rk( 31) * att
endif

rk( 32) =  dexp(-0.2590825451818744d+02                           &
 - ( -0.1800000000000000d+03 )/temp)
effko =  0.2500000000000000d-30* (temp / 3.d2)                    &
           **(- ( 0.1800000000000000d+01))
rapk =  0.2200000000000000d-10* (temp / 3.d2)                     &
            **(- ( 0.7000000000000000d+00))
rk( 33) = (effko * summ /                                         &
            ( 1.0d0 + effko * summ / rapk)) *                     &
            0.6000d+00** (1.0d0 / (1.0d0 +                        &
             (dlog10(effko * summ / rapk))**2))

if(option_photolysis.eq.2) then
 rk( 34)= 0.d0
elseif(option_photolysis.eq.1) then
if(azi.lt.0.d0)then
 stop
elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then
 rk( 34)=-0.1726633761595217d-06
 rk( 34)=-0.3763226238404795d-05+(azi- 0.00d+00) * rk( 34)
 rk( 34)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 34)
 rk( 34)= 0.2758968400000000d-01+(azi- 0.00d+00) * rk( 34)
elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then
 rk( 34)= 0.8277361284785716d-06
 rk( 34)=-0.8943127523190445d-05+(azi- 0.10d+02) * rk( 34)
 rk( 34)=-0.1270635376159524d-03+(azi- 0.10d+02) * rk( 34)
 rk( 34)= 0.2704069800000000d-01+(azi- 0.10d+02) * rk( 34)
elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then
 rk( 34)=-0.1363361137754772d-05
 rk( 34)= 0.1588895633116670d-04+(azi- 0.20d+02) * rk( 34)
 rk( 34)=-0.5760524953618990d-04+(azi- 0.20d+02) * rk( 34)
 rk( 34)= 0.2570348600000000d-01+(azi- 0.20d+02) * rk( 34)
elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then
 rk( 34)= 0.9462274225405240d-06
 rk( 34)=-0.2501187780147645d-04+(azi- 0.30d+02) * rk( 34)
 rk( 34)=-0.1488344642392872d-03+(azi- 0.30d+02) * rk( 34)
 rk( 34)= 0.2535296800000000d-01+(azi- 0.30d+02) * rk( 34)
elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then
 rk( 34)=-0.2183585524073524d-06
 rk( 34)= 0.3374944874739267d-05+(azi- 0.40d+02) * rk( 34)
 rk( 34)=-0.3652037935066590d-03+(azi- 0.40d+02) * rk( 34)
 rk( 34)= 0.2230966300000000d-01+(azi- 0.40d+02) * rk( 34)
elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then
 rk( 34)= 0.1378647870888786d-06
 rk( 34)=-0.3175811697481157d-05+(azi- 0.50d+02) * rk( 34)
 rk( 34)=-0.3632124617340762d-03+(azi- 0.50d+02) * rk( 34)
 rk( 34)= 0.1877676100000000d-01+(azi- 0.50d+02) * rk( 34)
elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then
 rk( 34)=-0.3865635959481289d-06
 rk( 34)= 0.9601319151840079d-06+(azi- 0.60d+02) * rk( 34)
 rk( 34)=-0.3853692595570321d-03+(azi- 0.60d+02) * rk( 34)
 rk( 34)= 0.1496492000000000d-01+(azi- 0.60d+02) * rk( 34)
elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then
 rk( 34)=-0.1916508555648358d-06
 rk( 34)=-0.1063677596325682d-04+(azi- 0.70d+02) * rk( 34)
 rk( 34)=-0.4821357000377863d-03+(azi- 0.70d+02) * rk( 34)
 rk( 34)= 0.1082067700000000d-01+(azi- 0.70d+02) * rk( 34)
elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then
 rk( 34)= 0.2540478756918318d-05
 rk( 34)=-0.1523639649680941d-04+(azi- 0.78d+02) * rk( 34)
 rk( 34)=-0.6891210797183578d-03+(azi- 0.78d+02) * rk( 34)
 rk( 34)= 0.6184712500000000d-02+(azi- 0.78d+02) * rk( 34)
elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then
 rk( 34)= 0.1651057353835306d-05
 rk( 34)= 0.4573509366928574d-04+(azi- 0.86d+02) * rk( 34)
 rk( 34)=-0.4451315023386027d-03+(azi- 0.86d+02) * rk( 34)
 rk( 34)= 0.9973396100000000d-03+(azi- 0.86d+02) * rk( 34)
elseif(azi.ge.90.d0)then
 rk( 34)= 0.5424277000000000d-04
endif
if(att.lt.0.99999) rk( 34) = rk( 34) * att
endif


return
end subroutine kinetic_2


!===============================================================================

!> fexchem_2
!> \brief Computation of the chemical production terms
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name              role
!------------------------------------------------------------------------------
!> \param[in]     ns                total number of chemical species
!> \param[in]     nr                total number of chemical reactions
!> \param[in]     y                 concentrations vector
!> \param[in]     rk                kinetic rates
!> \param[in]     zcsourc           source term
!> \param[in]     convers_factor    conversion factors
!> \param[out]    chem              chemical production terms for every species
!______________________________________________________________________________

subroutine fexchem_2(ns,nr,y,rk,zcsourc,convers_factor,chem)

implicit none

! Arguments

integer nr,ns
double precision rk(nr),y(ns),chem(ns),zcsourc(ns)
double precision convers_factor(ns)

! Local variables

integer i
double precision w(nr)
double precision conc(ns)

do i=1,ns
 chem(i)=0.d0
enddo

!     Conversion mg/kg to molecules/cm3.

do i = 1, ns
   conc(i) = y(i) * convers_factor(i)
enddo

!     Compute reaction rates.

call rates_2(ns,nr,rk,conc,w)

!     Chemical production terms.

chem( 19) = chem( 19) +  0.2000000000000000d+01 * w(  1)
chem( 20) = chem( 20) -  0.2000000000000000d+01 * w(  1)
chem( 16) = chem( 16) - w(  2)
chem( 19) = chem( 19) + w(  2)
chem( 20) = chem( 20) - w(  2)
chem( 17) = chem( 17) + w(  3)
chem( 19) = chem( 19) - w(  3)
chem( 20) = chem( 20) + w(  3)
chem( 16) = chem( 16) + w(  4)
chem( 17) = chem( 17) - w(  4)
chem( 16) = chem( 16) - w(  5)
chem( 17) = chem( 17) + w(  5)
chem(  2) = chem(  2) + w(  6)
chem( 16) = chem( 16) - w(  6)
chem(  2) = chem(  2) - w(  7)
chem( 17) = chem( 17) + w(  7)
chem(  2) = chem(  2) - w(  8)
chem( 15) = chem( 15) +  0.2000000000000000d+01 * w(  8)
chem(  8) = chem(  8) - w(  9)
chem( 15) = chem( 15) + w(  9)
chem( 20) = chem( 20) + w(  9)
chem(  3) = chem(  3) - w( 10)
chem( 15) = chem( 15) +  0.2000000000000000d+01 * w( 10)
chem(  4) = chem(  4) + w( 11)
chem( 10) = chem( 10) - w( 11)
chem( 14) = chem( 14) +  0.2000000000000000d+01 * w( 11)
chem(  4) = chem(  4) + w( 12)
chem( 10) = chem( 10) - w( 12)
chem(  1) = chem(  1) - w( 13)
chem(  4) = chem(  4) + w( 13)
chem( 10) = chem( 10) + w( 13)
chem( 13) = chem( 13) + w( 13)
chem( 14) = chem( 14) +  0.2000000000000000d+01 * w( 13)
chem( 11) = chem( 11) - w( 14)
chem( 12) = chem( 12) + w( 14)
chem( 19) = chem( 19) + w( 14)
chem( 11) = chem( 11) + w( 15)
chem( 12) = chem( 12) - w( 15)
chem( 19) = chem( 19) - w( 15)
chem( 10) = chem( 10) + w( 16)
chem( 12) = chem( 12) - w( 16)
chem( 13) = chem( 13) + w( 16)
chem( 14) = chem( 14) + w( 16)
chem( 19) = chem( 19) + w( 16)
chem( 20) = chem( 20) - w( 16)
chem( 14) = chem( 14) - w( 17)
chem( 15) = chem( 15) + w( 17)
chem( 19) = chem( 19) + w( 17)
chem( 20) = chem( 20) - w( 17)
chem(  8) = chem(  8) + w( 18)
chem( 15) = chem( 15) - w( 18)
chem( 20) = chem( 20) - w( 18)
chem(  9) = chem(  9) + w( 19)
chem( 15) = chem( 15) - w( 19)
chem( 19) = chem( 19) - w( 19)
chem(  5) = chem(  5) - w( 20)
chem(  6) = chem(  6) + w( 20)
chem( 14) = chem( 14) + w( 20)
chem( 15) = chem( 15) - w( 20)
chem( 16) = chem( 16) - w( 21)
chem( 18) = chem( 18) + w( 21)
chem( 19) = chem( 19) - w( 21)
chem( 18) = chem( 18) - w( 22)
chem( 19) = chem( 19) +  0.2000000000000000d+01 * w( 22)
chem( 20) = chem( 20) - w( 22)
chem( 18) = chem( 18) - w( 23)
chem( 20) = chem( 20) + w( 23)
chem(  7) = chem(  7) + w( 24)
chem( 18) = chem( 18) - w( 24)
chem( 19) = chem( 19) - w( 24)
chem(  7) = chem(  7) - w( 25)
chem( 18) = chem( 18) + w( 25)
chem( 19) = chem( 19) + w( 25)
chem(  7) = chem(  7) - w( 26)
chem(  9) = chem(  9) +  0.2000000000000000d+01 * w( 26)
chem( 13) = chem( 13) - w( 27)
chem( 19) = chem( 19) + w( 27)
chem( 20) = chem( 20) - w( 27)
chem( 13) = chem( 13) -  0.2000000000000000d+01 * w( 28)
chem(  3) = chem(  3) + w( 29)
chem( 14) = chem( 14) -  0.2000000000000000d+01 * w( 29)
chem(  3) = chem(  3) + w( 30)
chem( 14) = chem( 14) -  0.2000000000000000d+01 * w( 30)
chem( 17) = chem( 17) +  0.8900000000000000d+00 * w( 31)
chem( 18) = chem( 18) - w( 31)
chem( 19) = chem( 19) +  0.8900000000000000d+00 * w( 31)
chem( 20) = chem( 20) +  0.1100000000000000d+00 * w( 31)
chem( 17) = chem( 17) - w( 32)
chem( 19) = chem( 19) - w( 32)
chem( 20) = chem( 20) + w( 32)
chem( 17) = chem( 17) - w( 33)
chem( 18) = chem( 18) + w( 33)
chem( 19) = chem( 19) - w( 33)
chem(  7) = chem(  7) - w( 34)
chem( 18) = chem( 18) + w( 34)
chem( 19) = chem( 19) + w( 34)

!    Conversion molecules/cm3/sec to mg/kg/sec.

do i = 1, ns
   chem(i) = chem(i) / convers_factor(i)
enddo

!     Volumic source terms.

do i=1,ns
chem(i)=chem(i)+zcsourc(i)
enddo

return
end subroutine fexchem_2


!===============================================================================

!> jacdchemdc_2
!> \brief Computation of the Jacobian matrix for atmospheric chemistry
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name               role
!------------------------------------------------------------------------------
!> \param[in]     ns                 total number of chemical species
!> \param[in]     nr                 total number of chemical reactions
!> \param[in]     y                  concentrations vector
!> \param[in]     convers_factor     conversion factors of mug/m3 to
!>                                   molecules/cm3
!> \param[in]     convers_factor_jac conversion factors for the Jacobian matrix
!>                                   (Wmol(i)/Wmol(j))
!> \param[in]     rk                 kinetic rates
!> \param[out]    jacc               Jacobian matrix
!______________________________________________________________________________

subroutine jacdchemdc_2(ns,nr,y,convers_factor,                     &
                    convers_factor_jac,rk,jacc)

implicit none

! Arguments

integer nr,ns
double precision rk(nr),y(ns),jacc(ns,ns)
double precision convers_factor(ns)
double precision convers_factor_jac(ns,ns)

! Local variables

integer i,j
double precision dw(nr,ns)
double precision conc(ns)

do j=1,ns
 do i=1,ns
  jacc(i,j)=0.d0
 enddo
enddo

!     Conversion mg/kg to molecules/cm3.

do i = 1, ns
   conc(i) = y(i) * convers_factor(i)
enddo

call dratedc_2(ns,nr,rk,conc,dw)

jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw(  1, 20)
jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw(  1, 20)
jacc( 20, 20) = jacc( 20, 20)- 0.2000000000000000d+01*dw(  1, 20)
jacc( 20, 20) = jacc( 20, 20)- 0.2000000000000000d+01*dw(  1, 20)
jacc( 16, 16) = jacc( 16, 16) - dw(  2, 16)
jacc( 16, 20) = jacc( 16, 20) - dw(  2, 20)
jacc( 19, 16) = jacc( 19, 16) + dw(  2, 16)
jacc( 19, 20) = jacc( 19, 20) + dw(  2, 20)
jacc( 20, 16) = jacc( 20, 16) - dw(  2, 16)
jacc( 20, 20) = jacc( 20, 20) - dw(  2, 20)
jacc( 17, 19) = jacc( 17, 19) + dw(  3, 19)
jacc( 19, 19) = jacc( 19, 19) - dw(  3, 19)
jacc( 20, 19) = jacc( 20, 19) + dw(  3, 19)
jacc( 16, 17) = jacc( 16, 17) + dw(  4, 17)
jacc( 17, 17) = jacc( 17, 17) - dw(  4, 17)
jacc( 16, 16) = jacc( 16, 16) - dw(  5, 16)
jacc( 17, 16) = jacc( 17, 16) + dw(  5, 16)
jacc(  2, 16) = jacc(  2, 16) + dw(  6, 16)
jacc( 16, 16) = jacc( 16, 16) - dw(  6, 16)
jacc(  2,  2) = jacc(  2,  2) - dw(  7,  2)
jacc( 17,  2) = jacc( 17,  2) + dw(  7,  2)
jacc(  2,  2) = jacc(  2,  2) - dw(  8,  2)
jacc( 15,  2) = jacc( 15,  2)+ 0.2000000000000000d+01*dw(  8,  2)
jacc(  8,  8) = jacc(  8,  8) - dw(  9,  8)
jacc( 15,  8) = jacc( 15,  8) + dw(  9,  8)
jacc( 20,  8) = jacc( 20,  8) + dw(  9,  8)
jacc(  3,  3) = jacc(  3,  3) - dw( 10,  3)
jacc( 15,  3) = jacc( 15,  3)+ 0.2000000000000000d+01*dw( 10,  3)
jacc(  4, 10) = jacc(  4, 10) + dw( 11, 10)
jacc( 10, 10) = jacc( 10, 10) - dw( 11, 10)
jacc( 14, 10) = jacc( 14, 10)+ 0.2000000000000000d+01*dw( 11, 10)
jacc(  4, 10) = jacc(  4, 10) + dw( 12, 10)
jacc( 10, 10) = jacc( 10, 10) - dw( 12, 10)
jacc(  1,  1) = jacc(  1,  1) - dw( 13,  1)
jacc(  4,  1) = jacc(  4,  1) + dw( 13,  1)
jacc( 10,  1) = jacc( 10,  1) + dw( 13,  1)
jacc( 13,  1) = jacc( 13,  1) + dw( 13,  1)
jacc( 14,  1) = jacc( 14,  1)+ 0.2000000000000000d+01*dw( 13,  1)
jacc( 11, 11) = jacc( 11, 11) - dw( 14, 11)
jacc( 12, 11) = jacc( 12, 11) + dw( 14, 11)
jacc( 19, 11) = jacc( 19, 11) + dw( 14, 11)
jacc( 11, 12) = jacc( 11, 12) + dw( 15, 12)
jacc( 11, 19) = jacc( 11, 19) + dw( 15, 19)
jacc( 12, 12) = jacc( 12, 12) - dw( 15, 12)
jacc( 12, 19) = jacc( 12, 19) - dw( 15, 19)
jacc( 19, 12) = jacc( 19, 12) - dw( 15, 12)
jacc( 19, 19) = jacc( 19, 19) - dw( 15, 19)
jacc( 10, 12) = jacc( 10, 12) + dw( 16, 12)
jacc( 10, 20) = jacc( 10, 20) + dw( 16, 20)
jacc( 12, 12) = jacc( 12, 12) - dw( 16, 12)
jacc( 12, 20) = jacc( 12, 20) - dw( 16, 20)
jacc( 13, 12) = jacc( 13, 12) + dw( 16, 12)
jacc( 13, 20) = jacc( 13, 20) + dw( 16, 20)
jacc( 14, 12) = jacc( 14, 12) + dw( 16, 12)
jacc( 14, 20) = jacc( 14, 20) + dw( 16, 20)
jacc( 19, 12) = jacc( 19, 12) + dw( 16, 12)
jacc( 19, 20) = jacc( 19, 20) + dw( 16, 20)
jacc( 20, 12) = jacc( 20, 12) - dw( 16, 12)
jacc( 20, 20) = jacc( 20, 20) - dw( 16, 20)
jacc( 14, 14) = jacc( 14, 14) - dw( 17, 14)
jacc( 14, 20) = jacc( 14, 20) - dw( 17, 20)
jacc( 15, 14) = jacc( 15, 14) + dw( 17, 14)
jacc( 15, 20) = jacc( 15, 20) + dw( 17, 20)
jacc( 19, 14) = jacc( 19, 14) + dw( 17, 14)
jacc( 19, 20) = jacc( 19, 20) + dw( 17, 20)
jacc( 20, 14) = jacc( 20, 14) - dw( 17, 14)
jacc( 20, 20) = jacc( 20, 20) - dw( 17, 20)
jacc(  8, 15) = jacc(  8, 15) + dw( 18, 15)
jacc(  8, 20) = jacc(  8, 20) + dw( 18, 20)
jacc( 15, 15) = jacc( 15, 15) - dw( 18, 15)
jacc( 15, 20) = jacc( 15, 20) - dw( 18, 20)
jacc( 20, 15) = jacc( 20, 15) - dw( 18, 15)
jacc( 20, 20) = jacc( 20, 20) - dw( 18, 20)
jacc(  9, 19) = jacc(  9, 19) + dw( 19, 19)
jacc(  9, 15) = jacc(  9, 15) + dw( 19, 15)
jacc( 15, 19) = jacc( 15, 19) - dw( 19, 19)
jacc( 15, 15) = jacc( 15, 15) - dw( 19, 15)
jacc( 19, 19) = jacc( 19, 19) - dw( 19, 19)
jacc( 19, 15) = jacc( 19, 15) - dw( 19, 15)
jacc(  5,  5) = jacc(  5,  5) - dw( 20,  5)
jacc(  5, 15) = jacc(  5, 15) - dw( 20, 15)
jacc(  6,  5) = jacc(  6,  5) + dw( 20,  5)
jacc(  6, 15) = jacc(  6, 15) + dw( 20, 15)
jacc( 14,  5) = jacc( 14,  5) + dw( 20,  5)
jacc( 14, 15) = jacc( 14, 15) + dw( 20, 15)
jacc( 15,  5) = jacc( 15,  5) - dw( 20,  5)
jacc( 15, 15) = jacc( 15, 15) - dw( 20, 15)
jacc( 16, 19) = jacc( 16, 19) - dw( 21, 19)
jacc( 16, 16) = jacc( 16, 16) - dw( 21, 16)
jacc( 18, 19) = jacc( 18, 19) + dw( 21, 19)
jacc( 18, 16) = jacc( 18, 16) + dw( 21, 16)
jacc( 19, 19) = jacc( 19, 19) - dw( 21, 19)
jacc( 19, 16) = jacc( 19, 16) - dw( 21, 16)
jacc( 18, 18) = jacc( 18, 18) - dw( 22, 18)
jacc( 18, 20) = jacc( 18, 20) - dw( 22, 20)
jacc( 19, 18) = jacc( 19, 18)+ 0.2000000000000000d+01*dw( 22, 18)
jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 22, 20)
jacc( 20, 18) = jacc( 20, 18) - dw( 22, 18)
jacc( 20, 20) = jacc( 20, 20) - dw( 22, 20)
jacc( 18, 18) = jacc( 18, 18) - dw( 23, 18)
jacc( 18, 19) = jacc( 18, 19) - dw( 23, 19)
jacc( 20, 18) = jacc( 20, 18) + dw( 23, 18)
jacc( 20, 19) = jacc( 20, 19) + dw( 23, 19)
jacc(  7, 18) = jacc(  7, 18) + dw( 24, 18)
jacc(  7, 19) = jacc(  7, 19) + dw( 24, 19)
jacc( 18, 18) = jacc( 18, 18) - dw( 24, 18)
jacc( 18, 19) = jacc( 18, 19) - dw( 24, 19)
jacc( 19, 18) = jacc( 19, 18) - dw( 24, 18)
jacc( 19, 19) = jacc( 19, 19) - dw( 24, 19)
jacc(  7,  7) = jacc(  7,  7) - dw( 25,  7)
jacc( 18,  7) = jacc( 18,  7) + dw( 25,  7)
jacc( 19,  7) = jacc( 19,  7) + dw( 25,  7)
jacc(  7,  7) = jacc(  7,  7) - dw( 26,  7)
jacc(  9,  7) = jacc(  9,  7)+ 0.2000000000000000d+01*dw( 26,  7)
jacc( 13, 13) = jacc( 13, 13) - dw( 27, 13)
jacc( 13, 20) = jacc( 13, 20) - dw( 27, 20)
jacc( 19, 13) = jacc( 19, 13) + dw( 27, 13)
jacc( 19, 20) = jacc( 19, 20) + dw( 27, 20)
jacc( 20, 13) = jacc( 20, 13) - dw( 27, 13)
jacc( 20, 20) = jacc( 20, 20) - dw( 27, 20)
jacc( 13, 13) = jacc( 13, 13)- 0.2000000000000000d+01*dw( 28, 13)
jacc( 13, 13) = jacc( 13, 13)- 0.2000000000000000d+01*dw( 28, 13)
jacc(  3, 14) = jacc(  3, 14) + dw( 29, 14)
jacc(  3, 14) = jacc(  3, 14) + dw( 29, 14)
jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 29, 14)
jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 29, 14)
jacc(  3, 14) = jacc(  3, 14) + dw( 30, 14)
jacc(  3, 14) = jacc(  3, 14) + dw( 30, 14)
jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 30, 14)
jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 30, 14)
jacc( 17, 18) = jacc( 17, 18)+ 0.8900000000000000d+00*dw( 31, 18)
jacc( 18, 18) = jacc( 18, 18) - dw( 31, 18)
jacc( 19, 18) = jacc( 19, 18)+ 0.8900000000000000d+00*dw( 31, 18)
jacc( 20, 18) = jacc( 20, 18)+ 0.1100000000000000d+00*dw( 31, 18)
jacc( 17, 17) = jacc( 17, 17) - dw( 32, 17)
jacc( 17, 19) = jacc( 17, 19) - dw( 32, 19)
jacc( 19, 17) = jacc( 19, 17) - dw( 32, 17)
jacc( 19, 19) = jacc( 19, 19) - dw( 32, 19)
jacc( 20, 17) = jacc( 20, 17) + dw( 32, 17)
jacc( 20, 19) = jacc( 20, 19) + dw( 32, 19)
jacc( 17, 17) = jacc( 17, 17) - dw( 33, 17)
jacc( 17, 19) = jacc( 17, 19) - dw( 33, 19)
jacc( 18, 17) = jacc( 18, 17) + dw( 33, 17)
jacc( 18, 19) = jacc( 18, 19) + dw( 33, 19)
jacc( 19, 17) = jacc( 19, 17) - dw( 33, 17)
jacc( 19, 19) = jacc( 19, 19) - dw( 33, 19)
jacc(  7,  7) = jacc(  7,  7) - dw( 34,  7)
jacc( 18,  7) = jacc( 18,  7) + dw( 34,  7)
jacc( 19,  7) = jacc( 19,  7) + dw( 34,  7)

do j = 1, ns
   do i = 1, ns
      jacc(i,j) = jacc(i,j) * convers_factor_jac(i,j)
   enddo
enddo


return
end subroutine jacdchemdc_2


!===============================================================================

!> rates_2
!> \brief Computation of reaction rates
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name               role
!------------------------------------------------------------------------------
!> \param[in]     ns                 total number of chemical species
!> \param[in]     nr                 total number of chemical reactions
!> \param[in]     rk                 kinetic rates
!> \param[in]     y                  concentrations vector
!> \param[out]    w                  reaction rates
!______________________________________________________________________________

subroutine rates_2(ns,nr,rk,y,w)

implicit none

! Arguments

integer nr,ns
double precision rk(nr),y(ns)

! Local variables

double precision w(nr)

w(  1) =  rk(  1) * y( 20) * y( 20)
w(  2) =  rk(  2) * y( 16) * y( 20)
w(  3) =  rk(  3) * y( 19)
w(  4) =  rk(  4) * y( 17)
w(  5) =  rk(  5) * y( 16)
w(  6) =  rk(  6) * y( 16)
w(  7) =  rk(  7) * y(  2)
w(  8) =  rk(  8) * y(  2)
w(  9) =  rk(  9) * y(  8)
w( 10) =  rk( 10) * y(  3)
w( 11) =  rk( 11) * y( 10)
w( 12) =  rk( 12) * y( 10)
w( 13) =  rk( 13) * y(  1)
w( 14) =  rk( 14) * y( 11)
w( 15) =  rk( 15) * y( 12) * y( 19)
w( 16) =  rk( 16) * y( 12) * y( 20)
w( 17) =  rk( 17) * y( 14) * y( 20)
w( 18) =  rk( 18) * y( 15) * y( 20)
w( 19) =  rk( 19) * y( 19) * y( 15)
w( 20) =  rk( 20) * y(  5) * y( 15)
w( 21) =  rk( 21) * y( 19) * y( 16)
w( 22) =  rk( 22) * y( 18) * y( 20)
w( 23) =  rk( 23) * y( 18) * y( 19)
w( 24) =  rk( 24) * y( 18) * y( 19)
w( 25) =  rk( 25) * y(  7)
w( 26) =  rk( 26) * y(  7)
w( 27) =  rk( 27) * y( 13) * y( 20)
w( 28) =  rk( 28) * y( 13) * y( 13)
w( 29) =  rk( 29) * y( 14) * y( 14)
w( 30) =  rk( 30) * y( 14) * y( 14)
w( 31) =  rk( 31) * y( 18)
w( 32) =  rk( 32) * y( 17) * y( 19)
w( 33) =  rk( 33) * y( 17) * y( 19)
w( 34) =  rk( 34) * y(  7)

return
end subroutine rates_2


!===============================================================================

!> dratedc_2
!> \brief Computation of derivatives of reaction rates
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name               role
!------------------------------------------------------------------------------
!> \param[in]     nr                 total number of chemical reactions
!> \param[in]     ns                 total number of chemical species
!> \param[in]     rk                 kinetic rates
!> \param[in]     y                  concentrations vector
!> \param[out]    dw                 derivatives of reaction rates
!______________________________________________________________________________

subroutine dratedc_2(ns,nr,rk,y,dw)

implicit none

! Arguments

integer nr,ns
double precision rk(nr),y(ns)

! Local variables

double precision dw(nr,ns)

dw(  1, 20) =  rk(  1) * y( 20)
dw(  1, 20) =  rk(  1) * y( 20)
dw(  2, 16) =  rk(  2) * y( 20)
dw(  2, 20) =  rk(  2) * y( 16)
dw(  3, 19) =  rk(  3)
dw(  4, 17) =  rk(  4)
dw(  5, 16) =  rk(  5)
dw(  6, 16) =  rk(  6)
dw(  7,  2) =  rk(  7)
dw(  8,  2) =  rk(  8)
dw(  9,  8) =  rk(  9)
dw( 10,  3) =  rk( 10)
dw( 11, 10) =  rk( 11)
dw( 12, 10) =  rk( 12)
dw( 13,  1) =  rk( 13)
dw( 14, 11) =  rk( 14)
dw( 15, 12) =  rk( 15) * y( 19)
dw( 15, 19) =  rk( 15) * y( 12)
dw( 16, 12) =  rk( 16) * y( 20)
dw( 16, 20) =  rk( 16) * y( 12)
dw( 17, 14) =  rk( 17) * y( 20)
dw( 17, 20) =  rk( 17) * y( 14)
dw( 18, 15) =  rk( 18) * y( 20)
dw( 18, 20) =  rk( 18) * y( 15)
dw( 19, 19) =  rk( 19) * y( 15)
dw( 19, 15) =  rk( 19) * y( 19)
dw( 20,  5) =  rk( 20) * y( 15)
dw( 20, 15) =  rk( 20) * y(  5)
dw( 21, 19) =  rk( 21) * y( 16)
dw( 21, 16) =  rk( 21) * y( 19)
dw( 22, 18) =  rk( 22) * y( 20)
dw( 22, 20) =  rk( 22) * y( 18)
dw( 23, 18) =  rk( 23) * y( 19)
dw( 23, 19) =  rk( 23) * y( 18)
dw( 24, 18) =  rk( 24) * y( 19)
dw( 24, 19) =  rk( 24) * y( 18)
dw( 25,  7) =  rk( 25)
dw( 26,  7) =  rk( 26)
dw( 27, 13) =  rk( 27) * y( 20)
dw( 27, 20) =  rk( 27) * y( 13)
dw( 28, 13) =  rk( 28) * y( 13)
dw( 28, 13) =  rk( 28) * y( 13)
dw( 29, 14) =  rk( 29) * y( 14)
dw( 29, 14) =  rk( 29) * y( 14)
dw( 30, 14) =  rk( 30) * y( 14)
dw( 30, 14) =  rk( 30) * y( 14)
dw( 31, 18) =  rk( 31)
dw( 32, 17) =  rk( 32) * y( 19)
dw( 32, 19) =  rk( 32) * y( 17)
dw( 33, 17) =  rk( 33) * y( 19)
dw( 33, 19) =  rk( 33) * y( 17)
dw( 34,  7) =  rk( 34)

return
end subroutine dratedc_2


!===============================================================================

!> lu_decompose_2
!> \brief Computation of LU factorization of matrix m
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name               role
!------------------------------------------------------------------------------
!> \param[in]     ns                 matrix row number from the chemical species
!>                                   number
!> \param[in,out] m                  on entry, an invertible matrix.
!>                                   On exit, an LU factorization of m
!______________________________________________________________________________

subroutine lu_decompose_2 (ns,m)

implicit none

! Arguments

integer ns
double precision m(ns,ns)

! Local variables

double precision temp

!     Upper part.
m(2, 16) = m(2, 16) / m(2, 2)

!     Upper part.
m(3, 14) = m(3, 14) / m(3, 3)

!     Upper part.
m(4, 10) = m(4, 10) / m(4, 4)

!     Upper part.
m(5, 15) = m(5, 15) / m(5, 5)

!     Upper part.
temp = m(6, 5) * m(5, 15)
m(6, 15) = ( m(6, 15) - temp ) / m(6, 6)

!     Upper part.
m(7, 18) = m(7, 18) / m(7, 7)
!     Upper part.
m(7, 19) = m(7, 19) / m(7, 7)

!     Upper part.
m(8, 15) = m(8, 15) / m(8, 8)
!     Upper part.
m(8, 20) = m(8, 20) / m(8, 8)

!     Upper part.
m(9, 15) = m(9, 15) / m(9, 9)
!     Upper part.
temp = m(9, 7) * m(7, 18)
m(9, 18) = ( m(9, 18) - temp ) / m(9, 9)
!     Upper part.
temp = m(9, 7) * m(7, 19)
m(9, 19) = ( m(9, 19) - temp ) / m(9, 9)

!     Upper part.
m(10, 12) = m(10, 12) / m(10, 10)
!     Upper part.
m(10, 20) = m(10, 20) / m(10, 10)

!     Upper part.
m(11, 12) = m(11, 12) / m(11, 11)
!     Upper part.
m(11, 19) = m(11, 19) / m(11, 11)

!     Lower part.
temp = m(12, 11) * m(11, 12)
m(12, 12) = m(12, 12) - temp
!     Lower part.
temp = m(14, 10) * m(10, 12)
m(14, 12) = m(14, 12) - temp
!     Lower part.
temp = m(19, 11) * m(11, 12)
m(19, 12) = m(19, 12) - temp
!     Upper part.
temp = m(12, 11) * m(11, 19)
m(12, 19) = ( m(12, 19) - temp ) / m(12, 12)
!     Upper part.
m(12, 20) = m(12, 20) / m(12, 12)

!     Upper part.
temp = m(13, 12) * m(12, 19)
m(13, 19) = ( m(13, 19) - temp ) / m(13, 13)
!     Upper part.
temp = m(13, 12) * m(12, 20)
m(13, 20) = ( m(13, 20) - temp ) / m(13, 13)

!     Lower part.
temp = m(15, 3) * m(3, 14)
m(15, 14) = m(15, 14) - temp
!     Upper part.
temp = m(14, 5) * m(5, 15)
m(14, 15) = ( m(14, 15) - temp ) / m(14, 14)
!     Upper part.
temp = m(14, 12) * m(12, 19)
m(14, 19) = ( m(14, 19) - temp ) / m(14, 14)
!     Upper part.
temp = m(14, 10) * m(10, 20)
temp = temp + m(14, 12) * m(12, 20)
m(14, 20) = ( m(14, 20) - temp ) / m(14, 14)

!     Lower part.
temp = m(15, 5) * m(5, 15)
temp = temp + m(15, 8) * m(8, 15)
temp = temp + m(15, 14) * m(14, 15)
m(15, 15) = m(15, 15) - temp
!     Lower part.
temp = m(19, 14) * m(14, 15)
m(19, 15) = m(19, 15) - temp
!     Lower part.
temp = m(20, 8) * m(8, 15)
temp = temp + m(20, 14) * m(14, 15)
m(20, 15) = m(20, 15) - temp
!     Upper part.
temp = m(15, 2) * m(2, 16)
m(15, 16) = ( m(15, 16) - temp ) / m(15, 15)
!     Upper part.
temp = m(15, 14) * m(14, 19)
m(15, 19) = ( m(15, 19) - temp ) / m(15, 15)
!     Upper part.
temp = m(15, 8) * m(8, 20)
temp = temp + m(15, 14) * m(14, 20)
m(15, 20) = ( m(15, 20) - temp ) / m(15, 15)

!     Lower part.
temp = m(17, 2) * m(2, 16)
m(17, 16) = m(17, 16) - temp
!     Lower part.
temp = m(19, 15) * m(15, 16)
m(19, 16) = m(19, 16) - temp
!     Lower part.
temp = m(20, 15) * m(15, 16)
m(20, 16) = m(20, 16) - temp
!     Upper part.
m(16, 17) = m(16, 17) / m(16, 16)
!     Upper part.
m(16, 19) = m(16, 19) / m(16, 16)
!     Upper part.
m(16, 20) = m(16, 20) / m(16, 16)

!     Lower part.
temp = m(17, 16) * m(16, 17)
m(17, 17) = m(17, 17) - temp
!     Lower part.
temp = m(18, 16) * m(16, 17)
m(18, 17) = m(18, 17) - temp
!     Lower part.
temp = m(19, 16) * m(16, 17)
m(19, 17) = m(19, 17) - temp
!     Lower part.
temp = m(20, 16) * m(16, 17)
m(20, 17) = m(20, 17) - temp
!     Upper part.
m(17, 18) = m(17, 18) / m(17, 17)
!     Upper part.
temp = m(17, 16) * m(16, 19)
m(17, 19) = ( m(17, 19) - temp ) / m(17, 17)
!     Upper part.
temp = m(17, 16) * m(16, 20)
m(17, 20) = ( m(17, 20) - temp ) / m(17, 17)

!     Lower part.
temp = m(18, 7) * m(7, 18)
temp = temp + m(18, 17) * m(17, 18)
m(18, 18) = m(18, 18) - temp
!     Lower part.
temp = m(19, 7) * m(7, 18)
temp = temp + m(19, 17) * m(17, 18)
m(19, 18) = m(19, 18) - temp
!     Lower part.
temp = m(20, 17) * m(17, 18)
m(20, 18) = m(20, 18) - temp
!     Upper part.
temp = m(18, 7) * m(7, 19)
temp = temp + m(18, 16) * m(16, 19)
temp = temp + m(18, 17) * m(17, 19)
m(18, 19) = ( m(18, 19) - temp ) / m(18, 18)
!     Upper part.
temp = m(18, 16) * m(16, 20)
temp = temp + m(18, 17) * m(17, 20)
m(18, 20) = ( m(18, 20) - temp ) / m(18, 18)

!     Lower part.
temp = m(19, 7) * m(7, 19)
temp = temp + m(19, 11) * m(11, 19)
temp = temp + m(19, 12) * m(12, 19)
temp = temp + m(19, 13) * m(13, 19)
temp = temp + m(19, 14) * m(14, 19)
temp = temp + m(19, 15) * m(15, 19)
temp = temp + m(19, 16) * m(16, 19)
temp = temp + m(19, 17) * m(17, 19)
temp = temp + m(19, 18) * m(18, 19)
m(19, 19) = m(19, 19) - temp
!     Lower part.
temp = m(20, 12) * m(12, 19)
temp = temp + m(20, 13) * m(13, 19)
temp = temp + m(20, 14) * m(14, 19)
temp = temp + m(20, 15) * m(15, 19)
temp = temp + m(20, 16) * m(16, 19)
temp = temp + m(20, 17) * m(17, 19)
temp = temp + m(20, 18) * m(18, 19)
m(20, 19) = m(20, 19) - temp
!     Upper part.
temp = m(19, 12) * m(12, 20)
temp = temp + m(19, 13) * m(13, 20)
temp = temp + m(19, 14) * m(14, 20)
temp = temp + m(19, 15) * m(15, 20)
temp = temp + m(19, 16) * m(16, 20)
temp = temp + m(19, 17) * m(17, 20)
temp = temp + m(19, 18) * m(18, 20)
m(19, 20) = ( m(19, 20) - temp ) / m(19, 19)

!     Lower part.
temp = m(20, 8) * m(8, 20)
temp = temp + m(20, 12) * m(12, 20)
temp = temp + m(20, 13) * m(13, 20)
temp = temp + m(20, 14) * m(14, 20)
temp = temp + m(20, 15) * m(15, 20)
temp = temp + m(20, 16) * m(16, 20)
temp = temp + m(20, 17) * m(17, 20)
temp = temp + m(20, 18) * m(18, 20)
temp = temp + m(20, 19) * m(19, 20)
m(20, 20) = m(20, 20) - temp

return
end subroutine lu_decompose_2


!===============================================================================

!> lu_solve_2
!> \brief Resolution of MY=X where M is an LU factorization computed by
!>        lu_decompose_2
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------------
!   mode          name               role
!------------------------------------------------------------------------------
!> \param[in]     ns                 matrix row number from the chemical species number
!> \param[in]     m                  an LU factorization computed by lu_decompose_2
!> \param[in,out] x                  on entry, the right-hand side of the equation.
!                                    on exit, the solution of the equation
!______________________________________________________________________________

subroutine lu_solve_2 (ns, m, x)

implicit none

! Arguments

integer ns
double precision m(ns,ns)
double precision x(ns)

! Local variables

double precision temp

!     Forward substitution.

x(1) = x(1) / m(1, 1)

x(2) = x(2) / m(2, 2)

x(3) = x(3) / m(3, 3)

temp = m(4, 1) * x(1)
x(4) = ( x(4) - temp ) / m(4, 4)

x(5) = x(5) / m(5, 5)

temp = m(6, 5) * x(5)
x(6) = ( x(6) - temp ) / m(6, 6)

x(7) = x(7) / m(7, 7)

x(8) = x(8) / m(8, 8)

temp = m(9, 7) * x(7)
x(9) = ( x(9) - temp ) / m(9, 9)

temp = m(10, 1) * x(1)
x(10) = ( x(10) - temp ) / m(10, 10)

x(11) = x(11) / m(11, 11)

temp = m(12, 11) * x(11)
x(12) = ( x(12) - temp ) / m(12, 12)

temp = m(13, 1) * x(1)
temp = temp + m(13, 12) * x(12)
x(13) = ( x(13) - temp ) / m(13, 13)

temp = m(14, 1) * x(1)
temp = temp + m(14, 5) * x(5)
temp = temp + m(14, 10) * x(10)
temp = temp + m(14, 12) * x(12)
x(14) = ( x(14) - temp ) / m(14, 14)

temp = m(15, 2) * x(2)
temp = temp + m(15, 3) * x(3)
temp = temp + m(15, 5) * x(5)
temp = temp + m(15, 8) * x(8)
temp = temp + m(15, 14) * x(14)
x(15) = ( x(15) - temp ) / m(15, 15)

x(16) = x(16) / m(16, 16)

temp = m(17, 2) * x(2)
temp = temp + m(17, 16) * x(16)
x(17) = ( x(17) - temp ) / m(17, 17)

temp = m(18, 7) * x(7)
temp = temp + m(18, 16) * x(16)
temp = temp + m(18, 17) * x(17)
x(18) = ( x(18) - temp ) / m(18, 18)

temp = m(19, 7) * x(7)
temp = temp + m(19, 11) * x(11)
temp = temp + m(19, 12) * x(12)
temp = temp + m(19, 13) * x(13)
temp = temp + m(19, 14) * x(14)
temp = temp + m(19, 15) * x(15)
temp = temp + m(19, 16) * x(16)
temp = temp + m(19, 17) * x(17)
temp = temp + m(19, 18) * x(18)
x(19) = ( x(19) - temp ) / m(19, 19)

temp = m(20, 8) * x(8)
temp = temp + m(20, 12) * x(12)
temp = temp + m(20, 13) * x(13)
temp = temp + m(20, 14) * x(14)
temp = temp + m(20, 15) * x(15)
temp = temp + m(20, 16) * x(16)
temp = temp + m(20, 17) * x(17)
temp = temp + m(20, 18) * x(18)
temp = temp + m(20, 19) * x(19)
x(20) = ( x(20) - temp ) / m(20, 20)


!     Backward substitution.

temp = m(19, 20) * x(20)
x(19) = x(19) - temp

temp = m(18, 19) * x(19)
temp = temp + m(18, 20) * x(20)
x(18) = x(18) - temp

temp = m(17, 18) * x(18)
temp = temp + m(17, 19) * x(19)
temp = temp + m(17, 20) * x(20)
x(17) = x(17) - temp

temp = m(16, 17) * x(17)
temp = temp + m(16, 19) * x(19)
temp = temp + m(16, 20) * x(20)
x(16) = x(16) - temp

temp = m(15, 16) * x(16)
temp = temp + m(15, 19) * x(19)
temp = temp + m(15, 20) * x(20)
x(15) = x(15) - temp

temp = m(14, 15) * x(15)
temp = temp + m(14, 19) * x(19)
temp = temp + m(14, 20) * x(20)
x(14) = x(14) - temp

temp = m(13, 19) * x(19)
temp = temp + m(13, 20) * x(20)
x(13) = x(13) - temp

temp = m(12, 19) * x(19)
temp = temp + m(12, 20) * x(20)
x(12) = x(12) - temp

temp = m(11, 12) * x(12)
temp = temp + m(11, 19) * x(19)
x(11) = x(11) - temp

temp = m(10, 12) * x(12)
temp = temp + m(10, 20) * x(20)
x(10) = x(10) - temp

temp = m(9, 15) * x(15)
temp = temp + m(9, 18) * x(18)
temp = temp + m(9, 19) * x(19)
x(9) = x(9) - temp

temp = m(8, 15) * x(15)
temp = temp + m(8, 20) * x(20)
x(8) = x(8) - temp

temp = m(7, 18) * x(18)
temp = temp + m(7, 19) * x(19)
x(7) = x(7) - temp

temp = m(6, 15) * x(15)
x(6) = x(6) - temp

temp = m(5, 15) * x(15)
x(5) = x(5) - temp

temp = m(4, 10) * x(10)
x(4) = x(4) - temp

temp = m(3, 14) * x(14)
x(3) = x(3) - temp

temp = m(2, 16) * x(16)
x(2) = x(2) - temp

return
end subroutine lu_solve_2

